#===============================================================================
# Renshaw cell based on Traub, R. D. (1977) Repetitive firing of renshaw spinal
# interneurons. Biological Cybernetics, 27, 71--76.
# adapted to be independent from the surface area
#===============================================================================

from brian import *
from brian.library.random_processes import *
from brian.library.synapses import *

lengthD = 637 * um
lengthS = 106 * um
lengthIS = 35 * um

radiusD = 2.34 * um
radiusS = 2.34 * um
radiusIS = 1.5 * um

areaD = lengthD * 2 * radiusD * pi
areaS = lengthS * 2 * radiusS * pi
areaIS = lengthIS * 2 * radiusIS * pi

ratioS = 1
ratioD = areaD/areaS
ratioIS = areaIS/areaS

print ratioD
print ratioIS
Rin = 70 * ohm * cm

gleak = (1.0 / (3000.0 * ohm)) * cm ** -2

C = 1 * uF * cm ** -2

gnaS = 200.0 * msiemens * cm ** -2
gnaIS = 600.0 * msiemens * cm ** -2

gkS = 250 * msiemens * cm ** -2
gkIS = 200 * msiemens * cm ** -2

ENa = 115*mV
EK  = -5*mV


rcDtoS = (Rin * (lengthD / 2 + lengthS / 2)) / (radiusD ** 2 * pi)
gcDtoS = (1 / rcDtoS)/areaS

rcStoIS = (Rin * (lengthS / 2)) / (radiusS ** 2 * pi) + (Rin * (lengthIS / 2)) / (radiusIS ** 2 * pi)
gcStoIS = (1 / rcStoIS)/areaS


def g_inj(t):
    if t <= 4 * ms:
        return (0.00015 * sin((pi * t / ms) / 8) * msiemens)/areaS
    else:
        return (0.00015 * exp(-(t / ms - 4) / 20) * msiemens)/areaS
    return

def g_alpha(t):
    return t/ms*exp(-t/(10*ms))*5*msiemens*cm**-2

def g_biexp(t):
    return (exp(-t/(20*ms))-exp(-t/(1*ms)))*10*msiemens*cm**-2


eqs = '''alpha_mS = 1/ms*(0.8*(22.5-vm/mV)/(exp((22.5-vm/mV)/4)-1)) : 1/ms
        beta_mS = 1/ms*(0.7*(vm/mV-50)/(exp((vm/mV-50)/5)-1)) : 1/ms
        alpha_hS = 1/ms*(0.32*exp((35-vm/mV)/18)) : 1/ms
        beta_hS = 1/ms*(10/(exp((50-vm/mV)/5)+1)) : 1/ms
        alpha_nS = 1/ms*(0.03*(25-vm/mV)/(exp((25-vm/mV)/5)-1)) : 1/ms
        beta_nS = 1/ms*(0.5*exp((20-vm/mV)/40)) : 1/ms
        
        alpha_mIS = 1/ms*(0.8*(15-vIS/mV)/(exp((15-vIS/mV)/4)-1)) : 1/ms
        beta_mIS = 1/ms*(0.7*(vIS/mV-40)/(exp((vIS/mV-40)/5)-1)) : 1/ms
        alpha_hIS = 1/ms*(0.32*exp((40-vIS/mV)/18)) : 1/ms
        beta_hIS = 1/ms*(10/(exp((40-vIS/mV)/5)+1)) : 1/ms
        alpha_nIS = 1/ms*(0.03*(15-vIS/mV)/(exp((15-vIS/mV)/5)-1)) : 1/ms
        beta_nIS = 1/ms*(0.5*exp((10-vIS/mV)/40)) : 1/ms
        
        dmS/dt = alpha_mS*(1-mS)-beta_mS*mS : 1
        dhS/dt = alpha_hS*(1-hS)-beta_hS*hS : 1
        dnSo/dt = alpha_nS*(1-nSo)-beta_nS*nSo : 1
        
        dmIS/dt = alpha_mIS*(1-mIS)-beta_mIS*mIS : 1
        dhIS/dt = alpha_hIS*(1-hIS)-beta_hIS*hIS : 1
        dnIS/dt = alpha_nIS*(1-nIS)-beta_nIS*nIS : 1
        
        gepsp2 = g_alpha(t) : msiemens*cm**-2
        gepsp = epsp : msiemens*cm**-2
        
        i_inj : uA*cm**-2
        
        dvm/dt = -(i_inj+ gnaS*mS**3*hS*(vm-ENa) + gkS*nSo**4*(vm-EK) + gleak*vm + gcDtoS/ratioS*(vm-vD) + gcStoIS/ratioS*(vm-vIS))/C : mV
        dvIS/dt = -(gnaIS*mIS**3*hIS*(vIS-ENa) + gkIS*nIS**4*(vIS-EK) + gleak*vIS + gcStoIS/ratioIS*(vIS-vm))/C : mV
        dvD/dt = -(gepsp*(vD-60*mV)+gleak*vD+gcDtoS/ratioD*(vD-vm))/C : mV
        '''

tau1 = 1*ms
tau2 = 20*ms
eqs+=Equations('''dge/dt = -ge*invtau2 :siemens*metre**-2
                depsp/dt = (invpeak*ge-epsp)*invtau1 :msiemens*cm**-2
                dgi/dt = -gi*invtau2 :siemens*metre**-2
                dipsp/dt = (invpeak*gi-ipsp)*invtau1 :siemens*metre**-2  
                ''',invtau1=1/tau1,invtau2=1/tau2,invpeak = (tau2 / tau1) ** (tau1 / (tau2 - tau1)))

Gr = NeuronGroup(1, model=eqs,
threshold=EmpiricalThreshold(threshold= 40 * mV,refractory=5*ms),
implicit=True, clock=Clock(dt=0.01 * ms))

Gr.vm = 0 * mV
Gr.vD = 0 * mV
Gr.vIS = 0 * mV
Gr.mS = 0
Gr.hS = 1
Gr.nSo = 0
Gr.mIS = 0
Gr.hIS = 1
Gr.nIS = 0
Gr.ge= 8*msiemens*cm**-2#linspace(0.01,10,1000)*msiemens*cm**-2
Gr.i_inj=0*uA*cm**-2#-linspace(1,2000,1000)*uA*cm**-

tracevm = StateMonitor(Gr, 'vm', record=[0])
tracevD = StateMonitor(Gr, 'vD', record=[0])
tracevIS = StateMonitor(Gr, 'vIS', record=[0])
tracegepsp = StateMonitor(Gr, 'gepsp', record=[0])
traceepsp = StateMonitor(Gr, 'epsp', record=[0])

tracemS = StateMonitor(Gr, 'mS', record=[0])
tracehS = StateMonitor(Gr, 'hS', record=[0])
tracenSo = StateMonitor(Gr, 'nSo', record=[0])

spikesvm = SpikeMonitor(Gr)


figure(1)

run(120 * ms)
subplot(211)
tracevm.plot()

subplot(212)
tracegepsp.plot()
traceepsp.plot()
#tracevD.plot()
#subplot(313)
#tracevIS.plot()

figure(2)

subplot(311)
tracemS.plot()
subplot(312)
tracehS.plot()
subplot(313)
tracenSo.plot()

figure(3)
raster_plot(spikesvm)


#sptimes= spikesvm.getspiketimes()

#count =zeros(1000)

#for i in xrange(0,1000):
#    count[i] = sptimes[i].__len__()

#figure(4)
#plot(count,linspace(1,2000,1000))

show()
